Robustness of the avalanche dynamics in data packet transport on scale-free networks 
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We study the avalanche dynamics in the data packet transport on scale-free networks through a 
simple model. In the model, each vertex is assigned a capacity proportional to the load with the 
proportionality constant 1 + a. When the system is perturbed by a single vertex removal, the load 
of each vertex is redistributed, followed by subsequent failures of overloaded vertices. The avalanche 
size depends on the parameter a as well as which vertex triggers it. We find that there exists a 
critical value a c at which the avalanche size distribution follows a power law. The critical exponent 
associated with it appears to be robust as long as the degree exponent is between 2 and 3, and is 
close in value to that of the distribution of the diameter changes by single vertex removal. 

PACS numbers: 89.70.+C, 89.75.-k, 05.10.-a 



Avalanche dynamics, triggered by small initial pertur- 
bation, but spreading to other constituents successively, 
is one of intriguing problems in physics 0, 0, S El LA 
IE H> E3 EH TlTj . Such avalanche dynamics manifests 
itself in diverse forms such as cultural fads virus 
spreading 0, disease contagion 0, blackout in power 
transmissiongrids , data packet congestion in the 
Internet 0, 0, an( i so on - I n particular, the avalanche 
phenomena on complex networks are interesting, because 
they occur more frequently and their impact can be more 
severe than those occurring in the Euclidean space due to 
the close inter-connectivity among constituents in com- 
plex networks. 

To understand the intrinsic nature of the avalanche 
dynamics on complex networks, the sandpile model pro- 
posed by Bak, Tang, and Wiesenfeld has been studied on 
scale- free (SF) networks recently [l^. The SF network 
is the network whose degree distribution follows a power 
law, pd(k) ~ fc~ 7 . Since the sandpile model is a self- 
organized critical model, the avalanche size distribution 
follows a power law, p a (s) ~ s~ T , where s is the avalanche 
size. In the sandpile model, the exponent r depends on 
the degree exponent 7 of the embedded SF network as 
tbtw = 7/(7 — l) for 2 < 7 < 3 when the toppling 
threshold of each vertex is equal to its degree. However, 
when the toppling threshold is fixed as a constant, inde- 
pendent of degree, the exponent tmf = 3/2, being equal 
to the mean field value in the Euclidean space. Thus, 
it would be interesting to find an example of avalanche 
dynamics where the avalanche size distribution follows a 
power law with a nontrivial exponent, but different from 
the mean field value, and robust against variation of de- 
gree exponents. For this purpose, in this paper, we study 
the model proposed by Motter and Lai (ML) , designed 
to exploit the avalanche dynamics in the process of data 
packet transport on complex networks. 

In the ML model, each vertex is assigned a finite ca- 
pacity, given as 

Cj = (l + a)t?\ (1) 
where a is a control parameter and £j is the load of ver- 



tex j. The load of a given vertex is defined as the sum of 
the effective number of data packets passing through that 
vertex when every pair of vertices send and receive a unit 
data packet. The data packets are allowed to travel along 
the shortest pathways between a given pair of vertices 
and are divided evenly at each branching point |l3t 1 1 Jj ] . 
For SF networks, the load of each vertex is heterogeneous, 
and its distribution also follows a power law, pi{t) ~ t~ & ■ 
The superscript (0) in Eq. (1) indicates the load without 

any removal of vertices. The excess term a£j in Eq. (1) 
provides the ability to tolerate the additional burden and 
may describe the excess buffer at the routers in the In- 
ternet, for example. The basic assumption of the ML 
model is that the size of such an excess buffer is propor- 
tional to the activity at the vertex, the load £j. The 
control parameter a sets the global level of tolerance of 
the system. 

Next, we remove a vertex i intentionally, which we call 
triggering vertex. Then each pair of remaining vertices 
whose shortest pathway had passed through that trigger- 
ing vertex should find detours, resulting in rearrangement 
of the shortest pathways over the network, and the load 
at a remaining vertex j takes a new value, which is de- 
noted as £j . If the load £j exceeds its capacity Cj given 
by Eq. Q , then the vertex j would fail irreversibly. Other 
overloaded vertices also fail at the same time. These are 
the failures by the first shock, marked green with the 
symbol (o) in Fig. 1. After then, the shortest pathway 
configurations would rearrange again, and the overloaded 
vertices fail successively until no overloaded vertices re- 
main. The avalanche size s, is defined as the total num- 
ber of failed vertices throughout the whole process of 
the avalanche triggered by the vertex i. Note that in 
this model, failures do not necessarily proceed contigu- 
ously, that is, through the neighbors of vertices previ- 
ously failed, but spread over the entire system through 
nonlocal dynamics as shown in Fig. 1. For such nonlo- 
cal dynamics, the branching process formalism cannot be 
used to obtain the avalanche size distribution of the ML 
model. 
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FIG. 1: (Color online) Plot of the avalanche dynamics pat- 
tern at a c — 0.15 for a given small-size network. Cascading 
failures starting from the central vertex spread in a nonlocal 
way following the steps, Q ( re d), C 1 (green), and □ (blue). 



In the original work, ML measured the ratio G; = 
N!/N, where N and N[ are the numbers of vertices be- 
fore and after cascading failures, respectively, when the 
triggering vertex is i. Note that the avalanche size corre- 
sponds to Si = N—N?. ML found that Gi depends on the 
degree fcj of the triggering vertex i as well as the control 
parameter a. When a is large (small), the capacity of 
each vertex is large (small) , so that the number of failed 
vertices is small (large) and Gi is close to one (zero). 
Moreover, when the degree of the triggering vertex is 
large (small) , Gi is close to zero (one) , and the system is 
vulnerable (robust). Such numerical results suggest that 
there may occur a phase transition in the avalanche size. 
In this paper, we hnd numerically that indeed there exists 
a critical value a c at which the avalanche size distribution 
follows a power law, p a (s) ~ s~ T . We also study various 
features of the avalanche dynamics at the critical point. 

Let us first investigate the distribution of {si}, the 
avalanche size distribution p a (s). For large (small) a, the 
number of overloaded vertices is small (large), so that the 
avalanche size is finite (diverges) and the system may be 
considered as in a subcritical (supercritical) phase. We 
find that there exists a characteristic value a c between the 
two regimes, where the avalanche size distribution follows 
a power law, p a (s) ~ s~ r as shown in Fig. EI Numerical 
simulations are performed for the Barabasi- Albert (BA) 
model with different degree exponent values. We find 
that a c ~ 0.15, and r w 2.1(1), both of which are likely 
to be robust for different degree exponents 7 as long as 
2 < 7 < 3. While it is not manifest why such a robust 
behavior occurs in the avalanche size distribution, it is 
noteworthy to remind that other problems related to the 
shortest pathways such as the load distribution and the 
diameter change distribution are also likely to be robust. 
Thus the robustness of the avalanche size distribution 
may be caused by the notion of the shortest pathway. For 
7 > 3, however, p a (s) decays with exponent larger than 



FIG. 2: Plot of the avalanche size distribution for the BA 
model at a c = 0.15 with different 7 = 3.0 (0)> 2-6 (□), and 
2.2 (o). The mean degree is 4, and the system size is TV = 
10 4 . The dotted line has a slope —2.1, drawn for reference. 
The avalanche size distributions are obtained by deleting each 
vertex i in turn and measuring the respective avalanche size 
Si, then tabulating the histogram of Si, normalized by the 
number of triggering vertices TV. Inset: the avalanche size 
distribution(cumulative) under the same condition for 7 = 
3.0, but with a = 0.11 (top), 0.15 (middle), and 0.2 (bottom). 
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FIG. 3: Plot of the avalanche size distribution for the BA 
model with 7 = 3 by the first shock (□), compared with the 
avalanche size distribution including the entire process (O)- 
The slopes of dotted and dashed lines are —2.3 and —2.1, 
respectively, drawn for reference. 



r w 2.1 or exponentially depending on 7. The avalanche 
size distribution by the first shock behaves differently as 
Pa(sf) ~ s^ 2 ' 3 , which is shown in Fig. 3. We also check 



the avalanche size distribution for real world networks. 
For the yeast protein interaction network and the Inter- 
net, we obtain r w 2.3(1) and r « 1-8(1), respectively, 
as shown in Fig. ^ Note that the degree exponent of 
the yeast protein interaction network is 7 ~ 3.4 |16|. 
slightly larger than 3, thus the exponent r ~ 2.3 is some- 
what larger than 2.1(1) obtained in the BA model for 
2 < 7 < 3. 

The deviation of the exponent r for the Internet is 
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FIG. 4: Plot of the avalanche size distribution for the yeast 
protein interaction network (O) an d the Internet (□) at a c — 
0.15. The slopes of dotted and dashed lines are —1.8 and 
—2.3, respectively, drawn for reference. 



rooted from its different network structure from those of 
the protein interaction networks or the BA-type model 
networks: it was found through recent several studies 
that the Internet structure is effectively tree-like, while 
the protein interaction network and BA-type model net- 
works contain diverse connections [lj, • Accordingly, 
when a vertex on a branch of tree structure is removed, 
the giant cluster is divided into two or more components, 
and the giant cluster size shrinks apparently. Such a case 
occurs more often in the Internet than in other-type net- 
works, because the Internet is tree-like. Due to this fact, 
the avalanche size statistics of the Internet is different 
from that of other networks. On the other hand, while 
the rule of the ML model may not be relevant to the dy- 
namics in the protein interaction network or cellular net- 
works, cascading failure occurring in cellular networks is 
an important concept. For example, the protein interac- 
tion network provides the basic operational protocol in 
various signal transduction and functional pathways. In 
such a system, when a certain element (a protein or a 
substrate) fails or is removed (perturbed), others should 
take over its burden to survive the lack thereof, although 
the mechanism by which the cascade spreads could be 
different from that of the ML model. Studies in this di- 
rection have recently been carried out for the metabolic 
networks [H Ejj. 

Next, we examine the relationship of the mean 
avalanche size, denoted by (s) (fcdei), over different trig- 
gering vertices but with a given degree fcdei at a c in Fig. [3] 
We find that the quantity (s) (fcdei) increases with in- 
creasing kdeb However, there occur large fluctuations 
in (s)(fcdci), hi particular, for small fcdei- Note that if 
the orderings in the magnitude of (s)(fcdel) and fcdei are 
preserved, p a (s)ds = pd(k)dk and hence one has the re- 
lation (s)(fcdei) ~ ^dei _1 ^ T-1 ' ) - Indeed, Fig. exhibits 
such a behavior. To examine the fluctuations of (s) (fcdei) 
for given fcdei, we consider the distribution function q(x) 



FIG. 5: Plot of the mean avalanche size (s) versus the degree 
of the triggering vertex fcdei for the BA model (7 = 3) at a c . 
Data points (O) are averaged over different avalanche sizes 
triggered by the vertices with a given fcdei- The standard 
deviation of each data point is represented by a bar. The 
slope of the dotted line is the theoretical value 1.8, drawn 
for reference. Inset: Cumulative plot of the avalanche size 
distribution with the rescaled quantity x — s/k'^~ 1 ^ l ' T ~ 1 * > for 
fcdd = 3 (o), 5 (□), 10 (o), and 15 (A). 



of the avalanche sizes for given fcdei with a rescaled quan- 
tity, x — s/fefc 1 . Shown in the inset of Fig.[S]are 
the data of the cumulative distribution of (?cum(^) for dif- 
ferent fcdei, which collapse onto a single curve exhibiting 
a fat-tail behavior as q(x) ~ x~ 3 ' 2 for large x. 

Next, to study how much a given vertex with degree 
fc is vulnerable or robust under a random vertex failure, 
we count the number of failures rij of a vertex j out of 
N cascading events when each of N vertices acts as the 
triggering vertex. At this point, it is convenient to con- 
sider the random variables x? which take the value 1 if 
the vertex j topples due to the triggering vertex i and 
otherwise. In terms of x?, V\, x, 1 — and 5~\ x? = rij. 
Let /(fc) be the average of rij/N over the vertices with 
degree fc. Fig. shows the function /(fc) versus fc. It 
increases with increasing fc for small fc and exhibits a 
peak in the intermediate range of fc. For large fc, /(fc) 
is almost independent of fc. This result implies that the 
vertices with degree in the intermediate range are more 
vulnerable. Meanwhile, we note that the asymptotic 
value of /(fc) for large fc is ~ 0(1/N). This is because 
a vertex with large fc hardly fails through the cascading 
failure process triggered by other vertices, but fails triv- 
ially when itself acts as the triggering vertex. The peak 
at the intermediate value of fc in /(fc) can be understood 
in the following heuristic way. The load itself quantifies 
the level of traffic coming through the vertex, so we can 
expect that the higher the load is, the more excess traffic 
it would get by the breakdown of other vertices. On the 
other hand, since the excess capacity atp is assigned in 
a multiplicative way, the higher the load is, the larger 
room to accommodate the excess is, reducing the occa- 
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FIG. 6: Plot of the failure fraction / versus degree k at a c 
for the BA model (7 = 3) with N = 10 4 . Data points are 
logarithmically binned. Error bars represent the standard 
deviations for each bin. 




FIG. 7: Plot of the logarithm (with base 10) of the failure 
correlation function c(fe, fcdel) as functions of the degrees of 
the failed vertex k and of the triggering vertex fcdei- Data 
are logarithmically binned to reduce fluctuations. Simulation 
is performed for the BA model with 7 = 3, iV = 3000 and 
averaged over 10 configurations. 



sion to be overloaded. These two factors compete each 
other, generating a peak in the intermediate range of k 
in f(k). That implies the vertices with intermediate de- 
grees are more vulnerable. Such a behavior can also be 
seen in the information cascade model of Watts . 

This result is also reminiscent of the avalanche dy- 
namics of the sandpile model. The hubs, vertices with 
large degrees, play a role of the reservoir against fail- 
ures 01 ■ We also consider the failure correlation func- 
tion c(k, fcdci), defined as the average of x? with the con- 
straints kj — k and fcj = fcdeli ki denoting the degree of a 
vertex i. The darkest region in the bottom-right corner 



of Fig. 7 indicates that vertices with small degrees easily 
fail by the trigger of vertices with large degrees, whereas 
the reverse rarely happen, particularly for large kd c \, as 
manifested by the white region in the upper-left part of 
Fig. 7. 

It is interesting to notice that the avalanche size dis- 
tribution behaves similarly to the diameter change dis- 
tribution |20j |. Diameter is the average number of hops 
between every pair of vertices. Let S ^ be the diame- 
ter of a given network, where the superscript (0) means 
unperturbed network. When the network is perturbed 
by the removal of a vertex i, the diameter changes ac- 
cordingly, and the diameter of the remaining network 
is denoted as d^K Then the dimensionless quantity 
A, = (<fW — is measured for all i, and then 

its distribution function, composed of {A^}, behaves as 
Pdc(A) ~ A - '' for large A. The exponent £ was mea- 
sured to be £ w 2.2(1) for most artificial SF networks 
including the BA model, insensitive to the degree expo- 
nent 7 as long as 2 < 7 < 3, and £ f=s 2.3(1) for the yeast 
protein interaction network, but £ « 1.7(1) for the Inter- 
net. All the above values of the exponent £ are close to 
corresponding values of r for the avalanche size distribu- 
tion of the ML model. In addition, the exponents r and £ 
are also close in values to the load distribution exponent 
6 except for a few examples such as the Internet. Thus, 
it would be interesting to investigate the origin of such 
coincidences on a fundamental level. 

Finally, it is noteworthy that recently Zhao et al. [2lJ 
also studied the phase transition of the cascading failure 
for the ML model. They estimated the critical point to 
be a c w 0.1 by comparing the load distribution before 
and after the deletion of the hub. Their estimation is not 
inconsistent with our numerical estimation. However, the 
avalanche size distribution studied in this work provides 
a better criterion for the phase transition point. 

In conclusion, we have studied the avalanche dynamics 
in the model proposed by Motter and Lai, describing the 
data packet transport on SF networks. Depending on the 
model parameter a, which controls the magnitude of the 
capacity of each vertex, the pattern of avalanche dynam- 
ics can change. For small a, cascading failure spreads over 
the entire system, corresponding to supercritical behav- 
ior in avalanche dynamics. While, for large a, cascading 
failure is confined in a small region, and avalanche size 
follows a subcritical behavior. At the critical point a c , 
the avalanche size distribution follows a power law with 
exponent r. The exponent r seems to be robust for dif- 
ferent degree exponent 7 as long as 2 < 7 < 3, and is 
likely to be close to the exponent of the diameter change 
distribution. 
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